# created Wednesday, June 16, 2010 at 21:54
# author: Mikis Stasinopoulos
# based on functions in MASS 
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
drop1All <- function (object, 
                           scope,  
                           test = c("Chisq", "none"), 
                           k = 2, 
                           sorted = FALSE, 
                           trace = FALSE, 
                           ...) 
{
#-----------------------
  drop1.scope<-function (terms1, terms2) 
  {
    terms1 <- terms(terms1, "mu")
    f2 <- if (missing(terms2)) 
        numeric(0)
    else attr(terms(terms2, "mu"), "factor")
    factor.scope(attr(terms1, "factor"), list(drop = f2))$drop
  }
#-----------------------
    #if (!what %in% object$par) 
    #    stop(paste(what, "is not a parameter in the object", "\n"))
    tl <- attr(terms(object, "mu"), "term.labels")
    if (missing(scope)) 
        {scope <- drop1.scope(object)}
    else 
       {
        if (!is.character(scope)) 
            scope <- attr(terms(update.formula(formula(object, "mu"), scope), "mu"), 
                "term.labels")
        if (!all(match(scope, tl, FALSE))) 
            stop("scope is not a subset of term labels")
       }
    ns <- length(scope)
    ans <- matrix(nrow = ns + 1, ncol = 2, dimnames = list(c("<none>", 
        scope), c("df", "AIC")))
    ans[1, ] <- extractAIC(object, scale, k = k, ...)
    for (i in seq(ns)) 
      {
        tt <- scope[i]
        if (trace) 
            cat("trying -", tt, "\n")
        nfit <- update(object,  as.formula(paste("~ . -", tt)), what="All", evaluate = FALSE, trace=FALSE)
        
       # nfit <- eval.parent(nfit)
        nfit <- try(eval.parent(nfit), silent=TRUE)
              if (any(class(nfit)%in%"try-error"))
                 { 
                  cat("Model with term ", tt, "has failed \n")       
                  ans[i + 1, ] <- NA# extractAIC(object, scale, k = k, ...)          
                 }
            else ans[i + 1, ] <- extractAIC(nfit, scale, k = k, ...)
      }
    dfs <- ans[1, 1] - ans[, 1]
    dfs[1] <- NA
    aod <- data.frame(Df = dfs, AIC = ans[, 2])
    o <- if (sorted) 
        order(aod$AIC)
    else seq(along = aod$AIC)
    test <- match.arg(test)
    if (test == "Chisq") 
      {
        dev <- ans[, 2] - k * ans[, 1]
        dev <- dev - dev[1]
        dev[1] <- NA
        nas <- !is.na(dev)
        P <- dev
        P[nas] <- pchisq(dev[nas], dfs[nas], lower.tail = FALSE)
        aod[, c("LRT", "Pr(Chi)")] <- list(dev, P)
      }
    aod <- aod[o, ]
    head <- c("Single term deletions", "\nModel:", deparse(as.vector(formula(object))))
    class(aod) <- c("anova", "data.frame")
    attr(aod, "heading") <- head
    aod
}
#----------------------------------------------------------------------------------------
#  addterm
add1All   <- function (object, 
                           scope,
                           test = c("Chisq", "none" ), 
                              k = 2, 
                         sorted = FALSE, 
                          trace = FALSE, 
                           ...) 
{
#-----------------------------
 add.scope <- function (terms1, terms2 ) 
    {
    terms1 <- terms(terms1)
    terms2 <- terms(terms2)
    factor.scope(attr(terms1, "factor"), list(add = attr(terms2, "factor")))$add  
   }
 #---------------------------
    if (missing(scope) || is.null(scope)) 
        stop("no terms in scope")
    if (!is.character(scope))
        scope <- add.scope(object, terms(update.formula(formula(object, "mu"), scope)))

    if (!length(scope)) 
        stop("no terms in scope for adding to object")
    ns <- length(scope)
    ans <- matrix(nrow = ns + 1, ncol = 2, dimnames = list(c("<none>", 
        scope), c("df", "AIC")))
    ans[1, ] <- extractAIC(object, scale, k = k, ...)
    for (i in seq(ns)) 
      {
        tt <- scope[i]
        if (trace) 
            cat("trying +", tt, "\n")
        nfit <- update(object, as.formula(paste("~ . +", tt)), what="All", trace=FALSE, evaluate = FALSE)
        nfit <- try(eval.parent(nfit), silent=TRUE)
              if (any(class(nfit)%in%"try-error"))
                 { 
                  cat("Model with term ", tt, "has failed \n")       
                  ans[i + 1, ] <- NA# extractAIC(object, scale, k = k, ...)          
                 }
            else ans[i + 1, ] <- extractAIC(nfit, scale, k = k, ...)
      }
    dfs <- ans[, 1] - ans[1, 1]
    dfs[1] <- NA
    aod <- data.frame(Df = dfs, AIC = ans[, 2])
    o <- if (sorted) 
        order(aod$AIC)
    else seq(along = aod$AIC)
    test <- match.arg(test)
    if (test == "Chisq") 
      {
        dev <- ans[, 2] - k * ans[, 1]
        dev <- dev[1] - dev
        dev[1] <- NA
        nas <- !is.na(dev)
        P <- dev
        P[nas] <- pchisq(dev[nas], dfs[nas], lower.tail = FALSE)
        aod[, c("LRT", "Pr(Chi)")] <- list(dev, P)
      }
    aod <- aod[o, ]
    head <- c("Single term additions for", "\nModel:", deparse(as.vector(formula(object))))
    class(aod) <- c("anova", "data.frame")
    attr(aod, "heading") <- head
    aod
}
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
# Venable and Ripley step AIC function
#----------------------------------------------------------------------------------------
stepGAICAll.B <-function(object, 
                    scope, 
                direction = c("both", "backward", "forward"), 
                    trace = T, 
                     keep = NULL, 
                    steps = 1000,
                    scale = 0,
                        k = 2,
                       ...)                    
                    
{
#--------------------------------------------------------------------
    mydeviance <- function(x, ...) 
      {
        dev <- deviance(x)
        if (!is.null(dev)) 
            dev
        else extractAIC(x, k = 0)[2]
      }
#--------------------------------------------------------------------
    cut.string <- function(string) 
      {
        if (length(string) > 1) 
            string[-1] <- paste("\n", string[-1], sep = "")
        string
      }
#--------------------------------------------------------------------
    re.arrange <- function(keep) 
      {
        namr <- names(k1 <- keep[[1]])
        namc <- names(keep)
          nc <- length(keep)
          nr <- length(k1)
        array(unlist(keep, recursive = FALSE), c(nr, nc), list(namr, 
            namc))
      }
#--------------------------------------------------------------------
    step.results <- function(models, fit, object, usingCp = FALSE) #
      {
        change <- sapply(models, "[[", "change")
            rd <- sapply(models, "[[", "deviance")
            dd <- c(NA, abs(diff(rd)))
           rdf <- sapply(models, "[[", "df.resid")
           ddf <- c(NA, abs(diff(rdf)))
           AIC <- sapply(models, "[[", "AIC")
       heading <- c("Stepwise Model Path \nAnalysis of Deviance Table", 
            "\nInitial  Model:", deparse(as.vector(formula(object))), 
            "\nFinal   Model:", deparse(as.vector(formula(fit))), 
            "\n")
           aod <- if (usingCp) 
            data.frame(Step = change, Df = ddf, Deviance = dd, 
                "Resid. Df" = rdf, "Resid. Dev" = rd, Cp = AIC, 
                check.names = FALSE)
        else data.frame(Step = change, Df = ddf, Deviance = dd, 
            "Resid. Df" = rdf, "Resid. Dev" = rd, AIC = AIC, 
            check.names = FALSE)
        attr(aod, "heading") <- heading
        class(aod) <- c("Anova", "data.frame")
        fit$anova <- aod
        fit
      }
#---------------------------------------------------------------------- 
    Terms <- terms(object)
    #if ("mu")
    # {
           object$formula <- Terms 
      object$call$formula <- Terms
     #} 
   # else
   #  {
   #   object[[paste("sigma","formula",sep=".")]] <- Terms
   #   object[[paste("sigma","formula",sep=".")]][[2]]<-NULL
   #   if (paste("sigma", "formula", sep=".")%in%names(object$call)) 
   #     object$call[[paste(what,"formula",sep=".")]] <- formula(Terms)[-2]
   #  else ##this is when the sigma formula is not defined
   #     {object$call[[paste("sigma","formula",sep=".")]] <- formula(Terms)[-2]
   #     names(object$call)[length(names(object$call))]<-paste("sigma","formula",sep=".")
   #     }
   #  }     
           md <- missing(direction)
    direction <- match.arg(direction)
     backward <- direction == "both" | direction == "backward"
      forward <- direction == "both" | direction == "forward"
    if (missing(scope)) 
      {
        fdrop <- numeric(0)
         fadd <- attr(Terms, "factors")
        if (md) 
            forward <- FALSE
      }
    else 
      {
        if (is.list(scope)) 
          {
            fdrop <- if (!is.null(fdrop <- scope$lower)) 
                attr(terms(update.formula(formula(object, what="mu"), fdrop), what = "mu"), "factors")
            else numeric(0)
             fadd <- if (!is.null(fadd <- scope$upper)) 
                attr(terms(update.formula(formula(object, what="mu"), fadd), what = "mu"), "factors")
          }
        else 
          {
             fadd <- if (!is.null(fadd <- scope)) 
                attr(terms(update.formula(formula(object, what="mu"), scope), what = "mu" ), "factors")
            fdrop <- numeric(0)
          }
       }
    models <- vector("list", steps)
    if (!is.null(keep)) 
      keep.list <- vector("list", steps)
    if (is.list(object) && (nmm <- match("nobs", names(object), 
        0)) > 0) 
         n <- object[[nmm]]
    else n <- length(residuals(object))
     fit <- object
    bAIC <- extractAIC(fit, scale, k = k, ...)
     edf <- bAIC[1]
    bAIC <- bAIC[2]
    if (is.na(bAIC)) 
        stop("AIC is not defined for this model, so stepAIC cannot proceed")
      nm <- 1
   Terms <- terms(fit, "mu")
    if (trace)
        #cat("Distribution parameter: ", what, "\n") 
        cat("Start:  AIC=", format(round(bAIC, 2)), "\n", cut.string(deparse(as.vector(formula(fit, what="mu")))), 
            "\n\n")
    models[[nm]] <- list(deviance = mydeviance(fit), df.resid = n - 
        edf, change = "", AIC = bAIC)
    if (!is.null(keep)) 
        keep.list[[nm]] <- keep(fit, bAIC)
   usingCp <- FALSE
    while (steps > 0) 
     {
      steps <- steps - 1
        AIC <- bAIC
       ffac <- attr(Terms, "factors")
        if (!is.null(sp <- attr(Terms, "specials")) && !is.null(st <- sp$strata)) 
            ffac <- ffac[-st, ]
       scope <- factor.scope(ffac, list(add = fadd, drop = fdrop))
         aod <- NULL
      change <- NULL
        if (backward && length(scope$drop)) 
           {
            aod <- drop1All(fit, scope$drop, trace = max(0,trace - 1), k = k,test="none", ...)
             rn <- row.names(aod)
            row.names(aod) <- c(rn[1], paste("-", rn[-1], sep = " "))
            if (any(aod$Df == 0, na.rm = TRUE)) 
               {
                zdf <- aod$Df == 0 & !is.na(aod$Df)
                 nc <- match(c("Cp", "AIC"), names(aod))
                 nc <- nc[!is.na(nc)][1]
                 ch <- abs(aod[zdf, nc] - aod[1, nc]) > 0.01
                if (any(ch)) 
                  {
                  warning("0 df terms are changing AIC")
                  zdf <- zdf[!ch]
                  }
                if (length(zdf) > 0) 
                  change <- rev(rownames(aod)[zdf])[1]
               }
           }
        if (is.null(change)) 
           {
            if (forward && length(scope$add)) 
              {
                aodf <- add1All(fit, scope$add, trace = max(0, trace - 1), k = k,test="none", ...)
                  rn <- row.names(aodf)
                row.names(aodf) <- c(rn[1], paste("+", rn[-1], 
                  sep = " "))
                 aod <- if (is.null(aod)) 
                  aodf
                else rbind(aod, aodf[-1, , drop = FALSE])
              }
            attr(aod, "heading") <- NULL
            if (is.null(aod) || ncol(aod) == 0) 
                break
            nzdf <- if (!is.null(aod$Df)) 
                aod$Df != 0 | is.na(aod$Df)
             aod <- aod[nzdf, ]
            if (is.null(aod) || ncol(aod) == 0) 
                break
            nc <- match(c("Cp", "AIC"), names(aod))
            nc <- nc[!is.na(nc)][1]
             o <- order(aod[, nc])
            if (trace) 
                print(aod[o, ])
            if (o[1] == 1) 
                break
            change <- rownames(aod)[o[1]]
           }
        usingCp <- match("Cp", names(aod), 0) > 0
        fit <- update(fit, paste("~ .", change), evaluate = FALSE, what="All", trace = FALSE) #MS
        fit <- eval.parent(fit)
        if (is.list(fit) && (nmm <- match("nobs", names(fit), 
            0)) > 0) 
            nnew <- fit[[nmm]]
        else nnew <- length(residuals(fit))
        if (nnew != n) 
            stop("number of rows in use has changed: remove missing values?")
        Terms <- terms(fit, "mu")
         bAIC <- extractAIC(fit, scale, k = k, ...)
          edf <- bAIC[1]
         bAIC <- bAIC[2]
        if (trace) 
            cat("\nStep:  AIC=", format(round(bAIC, 2)), "\n", 
                cut.string(deparse(as.vector(formula(fit, "mu")))), 
                "\n\n")
        if (bAIC >= AIC + 1e-07) 
            break
        nm <- nm + 1
        models[[nm]] <- list(deviance = mydeviance(fit), df.resid = n - 
            edf, change = change, AIC = bAIC)
        if (!is.null(keep)) 
            keep.list[[nm]] <- keep(fit, bAIC)
    }
    if (!is.null(keep)) 
        fit$keep <- re.arrange(keep.list[seq(nm)])
    step.results(models = models[seq(nm)], fit, object, usingCp)
}
#----------------------------------------------------------------------------------------
